Ground state properties of fluxlines in a disordered environment 



OO 



Heiko Rieger 

HLRZ, Forschungszentrum Jiilich, 52425 Jiilich, Germany 
(June 30, 1997) 

A new numerical method to calculate exact ground states of multi-fluxline systems with quenched 
disorder is presented, which is based on the minimum cost flow algorithm from combinatorial opti- 
mization. We discuss several models that can be studied with this method including their specific 
implementations, physically relevant observables and results: 1) the A'^-hne model with A'' fluxlines 
(or directed polymers) in a d-dimensional environment with point and/or columnar disorder and 
hard or soft core repulsion; 2) the vortex glass model for a disordered superconductor in the strong 
screening limit and 3) the Sine-Gordon model with random pase shifts in the strong coupling limit. 
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Dirty type II superconductors in a magnetic field are 
the most intensively studied represerttatives of elastic 
manifolds in a disordered environmentcl. Their paradig- 
matic description consists in an ensemble of magnetic 
fluxlines (or vortexlines) interacting strongly with point 
and/or columnar defects and among themselves. This 
complicated multi-line situation is usually reduced to the 
study oLa single line, a directed polymer in a random 
mediumH, a problem that possesses deep connection, also 
to nonequilibrium fluctuations of moving interfacesEI. 

Here we are going to present a new numerical method 
(in the spirit of other recent applications of combinatorial 
optimization tools in the physics of disordered systemsa) 
by which the investigation of the full multi-line situa- 
tion becomes feasible. It will enable us to determine ex- 
act ground sates (i.e. minimum energy configurations) of 
theses systems in polynomial time. Since the low tem- 
perature physics of fluxlines in a random environment is 
dominated by disorder effects these ground state calcula- 
tions will enable us to make various statements about 
possibly glassy features, for instance the roughness of 
multi-line systems, the stiffness of vortex or gauge glass 
models and the displacement-displacement correlations 
in random phase models. 

To introduce the notation and to set the stage of the 
theoretical models we consider we start with a simple but 
non-trivial (and hence heavily discussedEI) example: the 
so called 1-line problem, which consists in determining 
the minimum energy configuration of a single (magnetic) 
fluxline or a directed polymer (for a Ill-lattice) in a dis- 
ordered environment. The lattice version of this model 
is given by the Hamiltonian (or energy function) 
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where is a sum over all bonds {ij) joining site i and 

j of a d-dimensional lattice, e.g. a rectangular {L'^~^ x H) 
lattice, with periodic boundary conditions (b.c.) in d — 1 
space direction and free b.c. in one direction. The bond 
energies eij > are quenched random variables that indi- 
cate how much energy it costs to put a segment of flux- 



line on a specific bond {ij). The fluxline configuration 
X {xij > 0), also called a flow, is given by specifying 
Xij = 1 for each bond i, which is occupied by the fluxline 



and X 



otherwise. For the configuration to form 
lines on each site of the lattice all incoming flow should 
balance the outgoing flow, i.e. the flow is divergence free 



V -x = , 



(2) 



where V- denotes the lattice divergence. Obviously the 
fluxline has to enter, and to leave, the system somewhere. 
We attach all sites of one free boundary to an extra site 
(via energetically neutral arcs, e = 0), which we call the 
source s, and the other side to another extra site, the 
target, t as indicated in fig. la. Now one can push one 
line through the system by inferring that s has a source 
strength of +1 and that t has a sink strength of —1, i.e. 



(V • x)^ = +N and (V • x)t 



-N 



(3) 



with N = 1. Thus, the 1-line problem consists in mini- 
mizing the energy (^ by finding a flow x in the network 
(the lattice plus the two extra sites s and t) fulfilling the 
constraints (||) and (^). 

The solution of this problem is equivalent to finding 
the shortest path from s to t, where distances between 
two lattice sites are identified with the energies Cij , which 
can either be done with Dijkstra's algorithm from com- 
binatorial optimizationLl or by equivalent methodfi. better 
known to physicists: the transfer matrix methodd. |-| 

Since this 1-line problem has been extensively studiedEI 
we directly proceed to its full generalization to N flux- 
lines, which has, to our knowledge, never been treated 
in the literature before. The reason is simple: Whereas 
two lines {N = 2) are still tractablecl, the transfer matrix 
method fails to work efficiently for an increasing number 
of lines since its complexity grows exponentially with N. 
Since it is the dense limit N = pL'^~^ with p of order one 
which is expected to contain new physics an algorithm 
that solves this problem in polynomial time as the one 
we are now going to present, is highly desirable. 

The A'^-line problem again consists in minimizing (|l|) 
in such a way that (|) and (|) are fulfilled, now with an 
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arbitrary value for N. Physically one has to take into 
account a repulsive interaction between the fluxlines, for 
instance a hard core repulsion, which can be modeled by 
inferring that Xij £ {0, 1}, i.e. that onlyj-a segment of one 
single fluxline can pass through an arda. It is also pos- 
sible to apply our method to a situation with soft core 
repulsion, which we discuss below. The problem is now 
formulated in such a way that it is identical to rainimum 
cost flow problem in combinatorial optimizatiorJj. 

Since the 1-line problem can be solved by finding a 
shortest path the intuitive idea to solve the iV-line prob- 
lem would be to search successively N shortest path, 
i.e. by adding one fluxline after the other to the system. 
However, adding a fluxline to an existing fluxline configu- 
ration might necessitate redirecting one or more fluxlines, 
as indicated in fig. 1. This at first sight formidable task is 
elegantly solved by the so called successive -shortest path 
algorithm for minimum cost flow problemscl. 




FIG. 1. Sketch of the residual network with shortest path 
before putting in th e first fluxline (top left) and after the up- 
date for one fluxline (top right). Bottom: Residual network 
after second iteration (left) and the actual optimal fluxline 
configuration (right). Note that this case is a non-trivial case 
when the first fluxline has to be deformed in order to achieve 
the optimal 2-line configuration. 

The first key ingredience is that one does not work 
with the original network but with the residual network 
corresponding to the actual fluxline configuration, which 
contains also the information about possibilities to send 
flow backwards (now with energy — Cy since one wins 
energy by reducing Xij), i.e. to modify the actual flow. 
Suppose that we put one fluxline along a shortest path 
P{s, t) from s to t, which means that we set Xij = 1 for 
all arcs on the path P(s, t). Then the residual network is 
obtained by reversing all arcs and inverting all energies 



along this path, indicating that here we cannot put any 
further flow in the forward direction (since we assume 
hard-core interaction, < 1), but can send flow 

backwards by reducing Xij on the forward arcs by one 
unit. This procedure is sketched in figure 1. 

The second key ingredience is the introduction of a so 
called potential (p that fulfills the relation 

VU) < Vii) + e«j- (4) 

for all arcs {ij) in the residual network, indicating how 
much energy ip{j) it would at least take to send one unit 
of flow from s to site j, IF it would cost an energy ip{i) 
to send it to site i. With the help of these potentials one 
defines the reduced costs 

cf^- = e,, + ^(z) - ^(j) > . (5) 

The last inequality, which follows from the properties of 
the potential ip <Q) actually ensures that there is no loop 
C in the current residual network (corresponding to a 
flow x) with negative total energy, since X](ijW£^ij — 
^ij^ implying that the flow x is optimaB. 

It is important to note that the inequality (jij) is rem- 
iniscent of a condition for shortest path distances d{i) 
from s to all sites i with respect to the energies e^j : they 
have to fulfill d{j) < d{i) + Cij. Thus, one uses these dis- 
tances d to construct the potential ip when putting one 
fluxline after the other into the network: 

We start with the empty network (zero fluxlines) x*^ = 
0, which is certainly an optimal flow for iV 0, and set 
(/? = 0, cfj = Cij. Next, let us suppose that we have an op- 
timal N — 1-line configuration corresponding to the fiow 
x^^^. The current potential is the reduced costs 

are c^~^ — Cij + tp^~^{i) — ^p^^^{j) and we consider the 
residual network G^~^ corresponding to the fiow x^~^ 
with the reduced costs c^^^ > 0. The iteration leading 
to an optimal iV-line configuration xfj is 

1. Determine a shortest path P{s,t) with respect to 
the reduced costs ci^~^from s to t in the residual 
network G^~^. 

2. For all site on P{s,t) let d{i) be the shortest path 
distance from s to i. For these update the poten- 
tials: cp^ii) = p'^-^i) + d{i) - d{t). 

3. To obtain xfj increase (decrease) by one unit the 
fiow variables x^j~^ on all forward (backward) arcs 
(ij) on the shortest path P{s,t). 

Note that due to the the fact that the numbers d{i) are 
shortest distances one has again cfj > 0, i.e. the flow x^ 
is indeed optimal. The complexity of this iteration is the 
same as that of Dijkstra's algorithm for finding shortest 
paths in a network, which is O(M^) in the worst case (M 
is the number of nodes in the network). We find, how- 
ever, for the cases we consider (d-dimensional lattices) it 
roughly scales linear in M = L"*. Thus, for N fluxlines 
the complexity of this algorithm is 0{NL^). 
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FIG. 2. Optimal fluxline configuration for one particular 
2d sample of linear size L = 100. The number of fluxlines, A^, 
is 1,2,4 (top), 8,16,32 (bottom). 
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FIG. 3. Scaling plot for the roughness (per line) in the 
two-dimensional multi-fluxline system, cf.(P). The data are 
averaged over 1000 different disorder configurations. 

In fig. 2 we show a number of optimal fluxline con- 
figurations with varying line density obtained with our 
algorithm for the A^-line problem on a L x L square lat- 
tice. We would like to emphasize the algorithm works for 
fluxlines in arbitrary dimension, even for general graphs. 
One observes that the roughness [Az«]av of the fluxlines 
decreases systematically with increasing line density from 
the single line limit {N =1) [Awjav ~ L'^ with the rough- 
ness exponent C = 2/3 to the dense limit {N = L in 
d = 2) with no roughness at all, [Au;]av = 0. For a fi- 
nite fluxline density p — N/L each fluxline is essentially 



free up to a length ^|| = ^^'^ with ^± = 1/p. 
expects for d = 2 the finite size scaling form 



[Aiujav L'^'w{Hp'^ ) with 



1/C 



Thus one 



(6) 



where H is the height and w is a scaling function with 
w{x) —^ const, for x ^ 0. In flg. 3 we show a correspond- 
ing scaling plot for thci-diata obtained with our algorithm 
for d = 2 (and H = Lf3. 



Before we proceed we would like to point out that the 
Hamiltonian (|l|) is general enough to describe various 
physically interesting situations. By an appropriate def- 
inition [Of the energies one can easily model columnar 
defectsQ, disorder induced melting of a fluxline lattice 
(e.g. the Abrikosov lattice) and depinning transitions. 

Next we consider soft core repulsion, which can be 
modeled by allowing a multiple occupancy of a bond 
(xij = 0, 1, 2, . . .) but punish high fluxline densities with 
an energy eij{xij) increasing faster than linear with the 
number of flux units Xij on the bond (ij). Thus the N- 
line problem with soft repulsion consists in minimizing 



i/(x) 



By {xi 



(7) 



under the constraints and (|3|) . The local energy func- 
tions Cij can be chosen arbitrarily for each bond (ij), 
however, they have to be convex as for instance Cij {xij ) — 
kij ■ x"j with n > 1 arbitrary. The energies have now 

°ij {-^ij ) ; 



to be replaced by the quantity eij{xij + 1) 
which is the energy needed to increase the flow Xij on 
arc (ij) by one unit. Since it depends on the current 
flow X the convexity of iij is needed to ensure that the 
reduced costs fulfill the inequality > (||) also after 
the flow modification. Whereas with hard core repulsion 
it was only possible to put N — L'^ fluxline into the sys- 
tem, the fluxline density can now arbitrarily high and an 
interplay between the repulsion and tlie disorder effects 
lead to a much richer phenomenologyllj. 

Up to now we considered situations in which the flux- 
lines are put into the system via an explicit external 
source. We now present a model in which fluxlines are 
generated inside the system: napaisly the vortex represen- 
tation of the gauge glass modelM strong screening, 
which is given by the HamiltonianLj 



— '^^{xij bij) . 



(8) 



Here ^ sum over all bonds of a simple cubic lat- 

tice {d — 3) with periodic b.c. in all directions and we 
do not have external source nodes subjected to condition 
(^. The Xij are the integer flow variables that have to 
fulfill the divergence free condition (H), and the bij are 
quenched random variables that are real numbers. They 
can be arbitrary, however in the gauge glass they fulfill a 
divergence free condition V • b = since they represent 
a magnetic field derived from a quenched random vector 
potential (b = V x A). Without the constraint (||) the 
optimal solution would simply be given by choosing Xij to 
be the closest integer to bij. This solution fulfills cfj > 
with (fi — 0, where the costs are chosen as for the convex 
fiow problem (Q). Since it violates the constraints (||) one 
has either excess or deficit flow entering or leaving indi- 
vidual sites, which one has to remove. Instead of sending 
flow from one particular source node s to a target t as 
in the fluxline problem, one now sends flow from excess 
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to deficit sites along shortest paths using the iteration 
described above. In this way one successively removes 
the violations of constraint (0) by keeping the reduced 
cost optimality cfj > all the time, which guarantees 
the optimality of the flow at the end of the iteration. 

The physically most interesting question in the context 
of the model Hamiltonian (^) concerns the existence of 
a superconducting glass phase. This can be-studied via 
domain wall renormalization group methodstJ, by which 
one determines the scaling behavior of low lying excita- 
tion AE on the length scale L. Such an excitation in 
the vortex representation (H) is a loop (closed fluxline) 
with an area proportional to L^, which can be realized by 
an extra fluxline (on the background of the true ground 
state) winding once around the 3d torus in one direction 
(note that we have periodic b.c. in all directions). De- 
tails of this procedure will be published elsewhere, here 
we give only the result: AE ~ with 6 = -0.95 ± 0.03 
From this we can draw two conclusions: a) since it is 
clearly negative there is no superconducting (vortex) 
glass phase at non-vanishing temperature, and b) the 
thermal correlation length diverges only with T ^ as 
Cth ~ T-" with an exponent l/\e\ - 1.05 ± 0.03, 
which is in agreement with Refs.E3. 

Finally, as our last application of our algorithm we dis- 
cuss the Sine-Gordon model with random phase shifts, 
which in two dimensions is a model for a fluxline arraytJ 

Hnp = - - cos (2^(«, - A)) , (9) 

(y) i 

where Ui are (real) displacements in a fluxline array. 
Pi € [0, 1[ (quenched) random phase shifts, {ij) nearest 
neighbor pairs of a c?-dimensional lattice, m > 1 is some 
number (usually it is m = 2, by which it the first term in 
(^ becomes an elastic energy) and A a coupling param- 
eter. In the limit of strong coupling A — > oo one enforces 
Ui = Pi + rii, where rii is an integer, and the Hamiltonian 
reads in these variables -ff^p = ~ ~ 

where 6ij = Oi — Oj. In d = 2 this model is equivalent j-to 
a model for a crystalline surface on a random substrateE3, 
for which the variables rii are interpreted as integer height 
variables and Oi as substrate heights between and 1. 

By introducing the variables Xij = rii— nj one observes 
that the form of the Hamiltonian is a special case of the 
convex cost functions considered above (^) for which the 
ground state can again computed with the algorithm we 
have presented. However, now one has to consider a sim- 
ply connected topology for the underlying network (e.g. 
not the torus) since we have to reconstruct the variables 
rii from their differences Xij. Thus we assume free or 
fixed b.c. for the lattice. In two dimensions d = 2 we 
studied the correlation function G{r) = [{rii — ?ij+r)^]av 
for various strengths of the non-linearity m, and found 
that G(r) increases stronger than log(r) with the dis- 
tance, possibly like log^(r), indicating a superrough low 
temperature phase, as for m = 2E3. 

To summarize we have presented various applications 



of a successive shortest path algorithm for disordered sys- 
tems containing many fluxlines or directed polymers with 
short range repulsion and derived a number of new re- 
sults. It would be of high interest to think about gen- 
eralization ta fluxlines with long range (like Coulomb) 
interaction.Ej 
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